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TIME DEPENDENCE OF COLLISION PROBABILITIES DURING 
SATELLITE CONJUNCTIONS 


Doyle T. Hall, Matthew D. Hejduk,t and Lauren C. Johnson? 


The NASA Conjunction Assessment Risk Analysis (CARA) team has recently 
implemented updated software to calculate the probability of collision (P.) for 
Earth-orbiting satellites. The algorithm can employ complex dynamical models 
for orbital motion, and account for the effects of non-linear trajectories as well as 
both position and velocity uncertainties. This “3D P,” method entails computing 
a 3-dimensional numerical integral for each estimated probability. Analysis indi- 
cates that the 3D method provides several new insights over the traditional “2D 
P.” method, even when approximating the orbital motion using the relatively sim- 
ple Keplerian two-body dynamical model. First, the formulation provides the 
means to estimate variations in the time derivative of the collision probability, or 
the probability rate, R.. For close-proximity satellites, such as those orbiting in 
formations or clusters, R- variations can show multiple peaks that repeat or blend 
with one another, providing insight into the ongoing temporal distribution of risk. 
For single, isolated conjunctions, R. analysis provides the means to identify and 
bound the times of peak collision risk. Additionally, analysis of multiple actual 
archived conjunctions demonstrates that the commonly used “2D P,” approxima- 
tion can occasionally provide inaccurate estimates. These include cases in which 
the 2D method yields negligibly small probabilities (e.g., P- < 10°!°), but the 3D 
estimates are sufficiently large to prompt increased monitoring or collision miti- 
gation (e.g., P- > 10°). Finally, the archive analysis indicates that a relatively 
efficient calculation can be used to identify which conjunctions will have negligi- 
bly small probabilities. This small-P. screening test can significantly speed the 
overall risk analysis computation for large numbers of conjunctions. 


INTRODUCTION 


As the population of artificial Earth-orbiting satellites grows, the likelihood of collisions be- 
tween satellites increases. Such collisions pose a direct threat to spacecraft occupied by humans, 
such as the International Space Station. They also threaten a much larger number of currently 
active robotic satellites, which represent significant investments by many international commercial, 
military and scientific organizations. Collisions also pose an indirect threat to future spacecraft, 
through the generation of fragmentation debris that can remain orbiting for extended periods. 
NASA’s CARA team has the responsibility to assess collision risks for a specific set of high-value 
satellites!, orbiting among the much larger population of cataloged objects’. Much of CARA’s 
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effort focuses on calculating the probability of collision during conjunction events, when two sat- 
ellites make a close approach to one another. Trajectories for cataloged satellites are not perfectly 
known; instead, they represent estimates based on tracking data gathered from a variety of sen- 
sors.*> Robust P. formulations must account for the uncertainties arising from this statistical orbital 
determination process. 


Previous Work 


Many authors have formulated analytical methods to estimate collision probabilities for Earth- 
orbiting satellites (see References 4-8, and references cited therein). Most of these employ some 
approximations that make the mathematical problem significantly more tractable, and ultimately 
allow P, values to be estimated using 2-dimensional numerical integration. Specifically, these “2D 
P.” methods approximate the curved relative satellite motion to be linear during the conjunction. 
They also assume that the duration of the conjunction is sufficiently short to allow the uncertainties 
on the satellite positions to be approximated as constant during the conjunction, and uncertainties 
on the satellite velocities to be neglected altogether. Approaches not using these approximations 
have also been discussed by several authors (e.g., References 6-10). Specifically, in 2012, V. T. 
Coppola developed an approach’ that estimates P. values using 3-dimensional numerical integra- 
tion. This “3D P.” method — which has been adopted by the CARA team because of its compre- 
hensive theoretical formulation — can use non-linear trajectories, and account for time-dependent 
uncertainties in both position and velocity. 


In 2015, K. Chan noted an issue with an assumption used in the 3D P- formulation,'! pointing 
out that Coppola’s integration of the flux of a time-dependent probability distribution function 
(PDF) through a hemi-spherical surface is not consistent with the basic tenets of probability theory. 
Because of this, there are limiting situations where the 3D P. method could potentially provide 
inaccurate results, which can be addressed by using an alternative analytical formulation based on 
the union of domains of random variables method.*""' In this alternative approach, P. values can be 
calculated using a 3D-volume integral within a carefully selected space of random variables. The 
integration volume comprises a tube swept out in this space by the collision sphere. In general this 
tube is curved. However, for most isolated satellite conjunctions, the tube curves gently relative to 
the size of the collision sphere or the PDF. In fact, the 2D P, approximation neglects this curvature 
altogether. Chan’s analysis'! indicates that Coppola’s method could provide inaccurate P. esti- 
mates in three general types of conjunctions: events where this tube crosses itself, events where the 
tube turns too sharply, or events where the PDF expands or contracts rapidly relative to the motion 
of the collision sphere. Preliminary analysis of more than 80,000 archived conjunctions reveals 
some relatively rare occurrences (<0.5%) of the first, tube-crossing event type, and none of the 
latter two types. As pointed out by both Coppola® and Chan,>!! such tube crossings can potentially 
occur in extended encounters between close-proximity satellites which produce multiple close ap- 
proach events. In these situations, P, estimates from Coppola’s method should be discarded in favor 
of Chan’s alternative method"! or a Monte Carlo method.* However, even in these situations, Cop- 
pola’s probability rate R. can still provide a useful indicator of how the collision risk varies in time. 
An example of such an extended encounter for two close-proximity satellites will be presented in 
this analysis. 


Summary of the 3D P, Implementation and Analysis 


This report explains the details of CARA’s recent implementation of Coppola’s 3D P, formula- 
tion, developed primarily as an incremental improvement to the current industry-standard 2D P. 
estimation method. This initial 3D P. implementation accounts for curved orbital motion as well 
as position and velocity uncertainties when estimating collision probabilities, and has been vali- 


dated using Monte Carlo simulations of several well-studied conjunctions.® This report also ex- 
plains how the 3D method can be used to estimate and analyze temporal variations in the collision 
probability rate,!* R., which is equal to the time derivative of P. for isolated conjunctions. Analysis 
of the probability rate provides the means to identify and bound the times of peak collision risk 
during isolated conjunction events much more precisely than the 2D P, estimation method. Anal- 
ysis of R, variations also provides the means to identify times of elevated risk for extended, blended, 
or repeating events, which occur for satellites persistently orbiting in close proximity (e.g., in clus- 
ters or formations). Such extended, blended, or repeating conjunction risks cannot be analyzed 
using 2D P, methods. Finally, analysis of archived conjunctions using the 3D P. method indicates 
that a relatively efficient pre-processing calculation can be used to identify which conjunctions will 
have negligibly small probabilities, providing the basis for a screening test that can significantly 
speed the risk analysis computation for large numbers of conjunctions. 


OVERVIEW OF COLLISION PROBABILITY THEORY 


This section provides a brief review of 2D and 3D P;, theory, focusing on how Coppola’s for- 
mulation® can be used to estimate variations in the collision probability rate, Re. 


Satellite Orbits, Uncertainties, and Probability Density Functions 


Earth-orbiting satellite trajectories are often modeled as “perturbed orbits” because the motion 
is affected by several perturbing forces including atmospheric drag, solar radiation pressure, solar 
and lunar gravity, and non-uniformities in Earth’s gravitational field.” Satellite tracking measure- 
ments enable estimation of the parameters defining such perturbed orbits.* For instance, the best- 
estimate trajectory for a satellite “s” may be described by specifying its Earth-centered inertial 
frame (ECI) position+velocity state vector,” x;(t-) = [1s(te) vs(t.)], at an initial time or epoch, f.. Best 
estimate states at other times, x,(t), can be propagated from this epoch state. (Note: in this discus- 


sion, such states will be considered to take the form of 6x1 column vectors.) 


Uncertainties on satellite states can be described using a 6-dimensional PDF. Such PDFs can 
also be estimated at a particular epoch using an orbit determination analysis of multiple tracking 
observations, and then propagated to a time (or set of times) of interest. ECI-state PDFs are often 
expressed using a “single-Gaussian” approximation. In this approximation, the PDF for satellite 
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s” can be written as a single multi-variate normal (MVN) function with dimensionality d = 6, as 


follows 
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where y denotes a possible state within the distribution, x; = x,(t) the center (or mean) of the distri- 
bution, P; = P,(f) the 6x6 state covariance matrix, and the superscript T indicates the transpose 
operator. The analysis presented here employs this single-Gaussian ECI-state PDF approximation 
for two reasons. First, this is the same approximation generally employed when using the 2D P. 
estimation method, and a large archive of actual historical conjunctions have been analyzed by the 
CARA team in this manner. So using this approximation in the earliest versions the 3D P, software 
exploits the single-Gaussian PDF descriptions readily available from archived data, and facilitates 
comparison of 2D and 3D results for a long record of historical conjunctions. Second, Coppola’s 
original formulation’® of the basic 3D P. equations also employs (but is not restricted to) such single- 
Gaussian PDFs. 


While the single-Gaussian approximation in Eq. (1) has been shown to represent a reasonable 
approximation for high-risk (P. > 10“) conjunctions,’ it can be regarded as only a first step towards 
amore general implementation of the 3D P, method. For instance, ECI-state PDFs can be approx- 
imated more accurately using a Gaussian mixture model (GMM)!*"!4:15 


G, 
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where i denotes the GMM summation index, w,, the weight of the i GMM component, x,,(f) the 
center of the i" component, P,,(t) the covariance matrix for the i" component, and G, the number 
of terms used in the GMM summation for satellite s. The 3D P. formulation can be extended to 
use such multi-component PDFs.° However, this generalization comes at the cost of significant 
notational complexity.!* For simplicity, the remainder of this discussion will only consider single 
Gaussian components. So the equations that follow can either be regarded as those appropriate for 
a single-Gaussian PDF approximation, or as isolated components within GMM summations (but 
with the associated summation indices suppressed for brevity). 


Finally, the size and shape of a six dimensional MVN distribution, N6 (y,x,P), is characterized 
by the 6x6 covariance matrix, P, which can be decomposed into three 3x3 sub-matrices® as follows 


A B’ 
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B C 
These sub-matrices A, B and C will be used later in this discussion. 


Relative Satellite States and Covariances 


The mean relative state for a pair of satellites (or for a pair of GMM components) can be defined 
as x(t) = [r(t) v()] = x;() — x,(t). Here, the subscripts p and s denote the primary and secondary 
satellites involved in a conjunction. The covariance for x(t) can be written P(t) = P,(t) + P,(d), 
assuming the two states x,(f) and x,(¢) are statistically independent.° The relative position covari- 
ance, A(t), as defined in Eq. (3), can be subjected to an Eigen-decomposition to estimate the orien- 
tation and size of the uncertainty ellipsoid for the relative position vector, r(#), of the two satellites.° 


Close Approach Events and Conjunctions 


Most conjunctions represent single, isolated close approach events between two satellites. 
These occur when the magnitude of the relative position vector, |r()|, reaches a minimum in time. 
The U.S. Air Force maintains orbital states for a large catalog of trackable satellites, enabling pre- 
dictions of such close approaches for CARA’s set of high-value primary satellites.’ Conjunction 
analysis occurs when any cataloged secondary satellite is predicted to make an incursion into a 
predefined screening volume centered on a primary.'''® For each conjunction, the CARA system 
processes an Orbital Conjunction Message (OCM) containing many quantities, including the time 
of close approach (t-2 or TCA). The OCM also contains data that can be readily converted to the 
ECI states for both the primary and secondary, Xp(fca) and X;(tca), respectively. These can be used 
to define the relative state at close approach, X(tca) = Xs(tca) — Xp(tca). Other OCM data can be readily 
converted to the ECI-state covariance matrices for the primary and secondary, P,(tca) and Ps(tca), 
respectively. Notably, the OCM does not contain the additional parameters required to describe 
the primary and secondary ECI PDFs as a multi-component GMM sum, as in Eq. (2). So using 
information from the OCM alone only readily enables the relative state PDF to be approximated 
using a single-Gaussian distribution, with ECI covariance P(tca) = Ps(tca) + Pp(tea). 


Approximations Used for 2D P, Estimation 


As mentioned previously, 2D P, estimation methods use some simplifying assumptions to ap- 
proximate P,.°° First, the curved relative orbital motion is approximated to be linear as follows 


r(t) <r(t,,)+[t—t, |v(t.) (4) 
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The 2D methods also assume that the duration of each conjunction is sufficiently short to allow two 
additional approximations: that the covariance of the relative position vector can be held constant 
during the conjunction, and that the uncertainties on the relative velocities can be neglected alto- 
gether, or 
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where 03;3 denotes a 3x3 matrix of zeros.° Using these approximations, the 2D methods can effi- 
ciently estimate a single P. value for each conjunction.*° 


Some conjunctions violate the assumptions behind these 2D P. approximations.”* For instance, 
not all conjunctions have a sufficiently short duration to justify their use. As discussed below, 
formulations to estimate conjunction durations*’ can be used to show that some events have ex- 
tremely long durations, comparable to the orbital periods of the primary or secondary satellites. 
Also, for satellite pairs persistently orbiting in close proximity, the concept of a single, isolated 
conjunction can break down altogether, because the events can repeat in a nearly periodic fashion, 
and become so extended that they blend into one another. 


Conjunction Bounds, Midpoints, and Durations 


This analysis estimates conjunction durations using a method also developed by Coppola.’ Alt- 
hough this formulation strictly applies to conjunctions analyzed using the 2D P. method, analysis 
indicates that it serves as a useful first-cut approximation for 3D conjunctions as well. The formu- 
lation’ denotes the beginning and ending bounding times of a conjunction as to(y) and t;(y), respec- 
tively, where y is a small number that defines a measure of achievable numerical precision, taken 
here to be y = 10°'°, roughly the resolution limit for computed double precision values. The con- 
junction duration spans these bounds: At = 1 — to. Similarly, the conjunction midpoint is defined 
here to bisect these bounds: tin = (t1 + To)/2. 


Other formulations to estimate conjunction durations*® have been demonstrated to be consistent 
with the method’ employed here. However, it should be reemphasized that these methods only 
strictly apply when using the short-encounter assumptions of the 2D P, method’ —i.e., the linearized 
dynamical model given in Eqs. (4)-(6). As will be discussed in more detail later, these conjunction 
durations and bounds do not always bracket the times of peak collision risk that occur when calcu- 
lating 3D P. values using more advanced dynamical models. In these cases, these first-cut con- 
junction bounds must be expanded appropriately. 


Formulation for 3D P, Estimation 


The 3D P, formulation (also referred to as the direct method") does not employ the restrictive 
2D-method assumptions discussed above, but instead can use relative ECI states and covariances 
from a dynamical model to estimate collision probabilities.° This dynamical model must provide 


the following set of time-dependent quantities: {x(), P(4)}, or alternatively {x(1), A, BY, C(O}. 
During the calculation, these quantities can be propagated from a reference epoch, or interpolated 
from a pre-tabulated ephemeris. The dynamical model could, conceivably, be very complex by 
including the effects of many perturbative forces. This approach would yield accurate results, but 
at the cost of considerable computation, some of which may not be required to achieve sufficient 
accuracy. Alternatively, the dynamical model could be very simple. For instance, nothing should 
prevent properly implemented 3D P. software from using the linearized dynamical model given in 
Eqs. (4)-(6). In this case, the 3D computation should always closely reproduce 2D P, values, as 
will be shown later. 


To demonstrate some of the differences that can arise between the 2D and 3D methods, this 
analysis also employs a slightly more complex dynamical model: Keplerian two-body orbital mo- 
tion? which describes elliptical trajectories within a uniform 1/r? gravitational field of a central 
body. This dynamical model still represents an approximation, because it neglects the effects of 
many important perturbative forces. However, over sufficiently short time intervals it does approx- 
imate the curved shape of the trajectories, and allows the full 6x6 covariances to be propagated in 
time.!” 


The expression for the 3D collision probability can be written as the sum of two components,° 
P. = Po + P;. Here, Po indicates the probability at an initial time, fo, and can often be neglected 
unless fo lies too close to a moment of peak collision risk. P; indicates the subsequent cumulative 
probability integrated over a finite time interval, fo < t < t9+ T. The CARA software explicitly 
calculates and includes Po when estimating 3D P. values. However, for simplicity, the remainder 
of this discussion will suppress Po, which allows the symbols P, and P; to be used interchangeably. 
With this simplification, P. could be expressed using a single equation containing a three-dimen- 
sional integral, as in Eq. (38b) of the original Coppola formulation.® However, following the for- 
mulation of DeMars er al '*, this analysis separates the expression into two equations in order to 
provide more insight into the time dependence of collision risk. The first equation contains the 
outermost integral over time 
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The second equation for the time-dependent probability rate,!* R.(f), contains the innermost two- 
dimensional integral, which spans the unit sphere 
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The integrand for this equation, denoted here using the symbol /, maps between points on the unit 
sphere and time, but also depends on the dynamical model quantities {x(4), P()}. This integrand 
comprises three factors, one of which is the square of the combined primary+secondary hard-body 
radii. The other two factors of J require several equations to describe fully.°!* These details are 
not repeated here. Instead, this discussion focuses on gaining a better understanding of the time 
dependence of collision risk provided by isolating the expression for the probability rate R.(f), and 
analyzing the associated temporal variations. For instance, as will be shown below, an isolated 
conjunction generally produces an R-(t) function with a single peak in time, which represents the 
moment of greatest collision risk. Satellites orbiting in close proximity, on the other hand, can 
produce multiple R, peaks. 


Monte Carlo P, Estimation 


Collision probabilities can also be estimated using Monte Carlo simulation methods.*!* These 
can be computationally intensive because they require repeatedly performing the following algo- 
rithmic steps: 1) sample the state PDFs of the primary and secondary satellites at some time; 2) 
propagate the sampled states over the conjunction time interval of interest, explicitly checking if 
the intervening distance ever becomes less than the combined hard-body radii; and 3) if so, register 
that a collision has occurred at the time of first contact between the two spheres defined by those 
radii. These steps need to be repeated until enough simulated collisions have been registered to 
provide sufficiently accurate statistical results. Step 2 requires the most computation, especially 
for complex dynamical models. However, for the linearized model in Eqs. (4) and (5), this step can 
be performed analytically, meaning that such simulations can be performed relatively quickly. 


In the sections below, calculations using the 3D method of Eqs. (7) and (8) will be compared to 
Monte Carlo simulations using both the linearized dynamical model and the more computationally 
intensive Keplerian two-body model. For these comparisons, the sampling performed in Monte 
Carlo step 1 above employs the single-Gaussian EC]-state distributions readily provided by the 
archived OCMs, as discussed above. Comparisons using more complex PDFs (e.g., multi-compo- 
nent GMMs) and/or more advanced dynamical models (e.g., high-fidelity state+covariance propa- 
gators) will be deferred for future analyses. 


CARA IMPLEMENTATION OF 3D Pc ESTIMATION 


This section outlines the details of how CARA’s 3D P- software calculates the integrals in Eqs. 
(7) and (8) numerically, as well as how the Keplerian two-body dynamical model has been imple- 
mented as a first step towards a more general and accurate approach. 


Numerical Integration over Time 


The 1D integral over time in Eq. (7) can be calculated using numerical techniques such as Simp- 
son’s rule integration® or trapezoidal integration. Both of these schemes require that the integration 
limits and step sizes be adjusted to achieve a desired numerical accuracy. As mentioned previously, 
the conjunction bounds (to, t1) only provide a first-cut approximation for the integration limits 
when using the Keplerian two-body dynamical model. For instance, the archive study indicates 
that, even for isolated conjunctions, (to, t1) do not always necessarily bracket the time when R.(f) 
peaks. Also, for archived conjunctions involving close-proximity satellites, (to, t1) can fail to 
bracket one or more of the repeating peaks in R.(t). Ideally, these first-cut integration limits should 
be adjusted as required during the computation to account for these effects, dynamically refined in 
order to find and sufficiently bracket all of the numerically significant peaks in R.(f). Similarly, 
the integration step size should also be dynamically adjusted. The CARA team plans to implement 
such adaptive integration schemes in future versions of the software. Currently, however, the con- 
junction bounds are simply expanded by a user-specified factor, E, resulting in the following lower 
and upper integration limits 
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where At and 1,, denote the conjunction duration and midpoint defined previously. When using 
the linearized dynamical model from Eqs. (4)-(6), an expansion factor of E = 1 should always 
suffice for these integration limits.’ Preliminary analysis indicates that for the Keplerian two-body 
dynamical model, E factors of 2 to 5, along with an integration step size of h = At/100, often suffice 
for single, isolated conjunctions. More conservatively, the expansion can be increased to E = 10, 
but often at the cost of more computation with no significant increase in integration accuracy. 
However, this simplified expansion approach does not address cases where there are multiple peaks 
in R(t), where, as mentioned previously, the concept of a single conjunction with well-defined 
temporal limits can break down altogether. Also, some cases have been found where the integration 
step size of h x At/100 seems to be too coarse, causing Simpson’s rule integrations to show unstable 
convergence behavior; for this reason, the current version of CARA’s software uses trapezoidal 
integration, which is less susceptible to such instabilities. 


Numerical Integration over the Unit Sphere 


The 2D unit sphere integral in Eq. (8) can be evaluated efficiently using a scheme called Lebe- 
dev quadrature.*!”'8 This method computes the 2D integral using a single weighted summation 
over a set of pre-defined points on the unit sphere. Quadrature points and weights have been de- 
fined'® for many different spherical harmonic or algebraic orders, up to a maximum of 131. In- 
creasing the algebraic order provides better integration accuracy at the cost of increased computa- 
tion. The CARA team plans eventually to implement a dynamic, self-adjusting integration scheme 
for this unit sphere integration. However, the current version of CARA’s software simply uses the 
maximum available Lebedev algebraic order of 131, which entails a weighted sum over 5810 
points. 


Covariance Propagation for Keplerian Two-Body Motion 


When using the Keplerian two-body dynamical model, CARA’s 3D P-. software currently propa- 
gates ECI-frame covariances using analytically-derived state transition matrices.'!’ More specifi- 
cally, the software employs the equation P(t) = ® P(t.a) ®’, where ® = (t, f-a) indicates the state- 
deviation transition matrix for two-body motion published by Shepperd.'’ This method of covari- 
ance propagation could be improved, because it can fail to produce realistic estimates for P(f) over 
long propagation periods.'”!?”° For instance, propagating covariances expressed in equinoctial or- 
bital elements may be a better approach’? especially for more complex dynamical models.”° How- 
ever, one of the goals of the current analysis is to compare the 2D and 3D P. methods, and investi- 
gate differences that can potentially arise by introducing nonlinear motion and time-dependent 6x6 
covariances. This goal can be accomplished using Shepperd’s state transition matrix formulation. 
However, once again, the CARA team regards this as only a first step towards a more general and 
robust approach. 


VALIDATION USING WELL-STUDIED CONJUNCTIONS 


In 2009, Alfano published a discussion of twelve conjunctions analyzed in detail using Monte 
Carlo and other P, estimation methods.* This includes one case where the 2D and Monte Carlo P- 
methods match one another very well, but also several cases where the 2D P, method fails to re- 
produce the Monte Carlo calculations. These well-studied conjunctions provide the means to val- 
idate the implementation of CARA’s new 3D P, software. Specifically, Alfano’s Monte Carlo 
analyses for these cases provide benchmark P, values that should be closely reproduced by the 3D 
P. software, if properly implemented. The discussion below focuses on two of these benchmark 
cases, selected to illustrate that the newly-implemented 3D P. software does indeed reproduce 
Alfano’s results (as it also does for the other benchmark cases not discussed or shown here). 


Validation Using Alfano’s Test Case 3 


Figure 1 shows the results of an analysis of “linear” test case 3, demonstrated by Alfano to be 
a conjunction where the 2D and Monte Carlo P, methods match one another.® The left panel 
plots the probability rate as a function of time, and the right panel plots the associated cumulative 
probability. Legends above each panel describe the plotted quantities and calculated P, values. 
The black points in each panel show the results of CARA’s Monte Carlo simulation using the lin- 
earized dynamical model from Eqs. (4)—(6), with a sample size of 30 million. The green curves 
show the output from CARA’s 3D P, software, also using the linearized dynamical model and an 
expansion factor of E = 1. The vertical cyan lines in the left panel show the conjunction bounds 
(To, T1), Which bracket the peak in probability rate (and which serve as the temporal integration 
limits in this case because EF = 1). The horizontal magenta line in the right panel shows the calcu- 
lated 2D P, value (labeled in the legend as “Foster P.’”’, in reference to the first author of an early 
formulation of the 2D method*). As discussed previously, when the 3D P. software is driven us- 
ing these linearized inputs, it should always closely reproduce the 2D P. value; the convergence 
of the green curve and the magenta line in the right panel demonstrates this explicitly. For this 
case, Alfano’s Monte Carlo simulation using 30 million samples yields P. = 0.10034, essentially 
identical to CARA’s 3D and Monte Carlo P, estimates cited in the figure legend. 
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Figure 1. Linearized-mode validation plots for Alfano’s conjunction test case 3. 


Step-by-Step Relaxation of the 2D P. Assumptions 


The assumptions employed in the 2D P, approximation can be relaxed in a step-by-step fashion, 
which provides additional insight when validating cases for which the 2D and 3D methods disagree. 
For the first step, the 3D P, software uses the linearized dynamical model from Eqs. (4)-(6), where 
it should always reproduce the 2D P, values, as demonstrated in the previous section. This first 
step can be described succinctly using the identifier “Linear motion, A=A(TCA), B=C=0,” notation 
which will appear on the legends of several figures presented later. For the second relaxation step, 
the linear motion can be replaced with Keplerian two-body motion, but the approximations of Eq. 
(6) left intact, described succinctly as “Kep2Body, A=A(TCA), B=C=0.” The third step introduces 
the time dependence of the relative position covariance, or “Kep2Body, A=A(t), B=C=0.”. Finally, 
the fourth step introduces the full 6x6 time-dependent covariance matrix, or “Kep2Body, P = P(/).” 
These relaxation steps define four operating modes of CARA’s current 3D P, software, identified 


here as “Coppola 1” through “Coppola 4,” respectively. Among these, the fourth mode “Coppola 
4” provides the current best estimate 3D P, values, to be compared to Alfano’s Monte Carlo bench- 
mark P, values for this validation study. 
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Figure 2. Step-by-step relaxation plots for Alfano’s conjunction test case 3. 


Figure 2 shows such a step-by-step relaxation analysis for Alfano’s linear test case 3, the same 
case shown in Figure 1. Again the two panels plot R- (left) and cumulative P. (right) against time. 
Curves for all four relaxation steps represent output from the 3D P. software, with the first step 
using an expansion factor of E = 1 and the other three E = 2. Legends above each panel show the 
plotted quantities and calculated P, values. The curves for all four of these relaxation steps plot 
directly on top of one another, and match the Foster 2D result, as they should for this benchmark 
case where the 2D and 3D methods are known to agree.® 


Validation Using Alfano’s Test Case 10 


Figure 3 shows a step-by-step relaxation analysis of Alfano’s “nonlinear” test case 10, using the 
same format as Figure 2. Again, curves for all four steps represent output from the 3D P. software, 
with the first step using an expansion factor of EF = 1, but the other three requiring E' ~ 5 for satis- 
factory integration convergence for this long duration conjunction. These plots clearly show that 
the four relaxation steps each produce distinctly different results. This is not surprising, because 
Alfano’s analysis* demonstrates that for this case the 2D and Monte Carlo P, values disagree, in- 
dicating that one or more of the assumptions used in the 2D method has been violated. Alfano’s 
Monte Carlo simulation using 1 billion samples yields a benchmark P, = 0.36300, differing by only 
0.3% from the best estimate 3D P, value of 0.36406. Both of these results differ significantly from 
the 2D estimate of P. = 0.29016. 
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Figure 3. Step-by-step relaxation plots for Alfano’s conjunction test case 10. 


Figure 4 shows this same test case 10 analyzed using CARA’s Monte Carlo simulation software 
with 3 million samples (orange diamonds) compared to the best estimate results from the 3D P. 
software (dotted black curve). In this case, the temporal integration limits have been adjusted to 
match those used in Alfano’s original analysis.’ Clearly, CARA’s Monte Carlo simulation closely 
reproduces the 3D P-. software results, and both match Alfano’s benchmark P- value to within 0.3%. 
The left panel of Figure 4 clearly indicates that this long duration conjunction comprises two R-(f) 
peaks which occur sufficiently close in time that they blend together. Inspection of the cumulative 
P. plotted as the dotted black curve in the right panel indicates that the 3D P, calculation reproduces 
Alfano’s corresponding benchmark curve (see Figure A10 of Reference 8). 
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Figure 4. Monte Carlo and 3D method analysis for Alfano’s conjunction test case 10. 
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ANALYSIS OF ARCHIVED CONJUNCTION EVENTS 


This section describes the preliminary analysis of some actual historical conjunctions previously 
analyzed by the CARA team. Specifically, the analysis employs a set of 80,827 archived OCMs, 
representing conjunctions that occurred between 2016-04-01 and 2016-06-19. Of these, 80,453 
are evaluated to be “high quality” events, characterized by positive-definite covariances A(tea) and 


PG): 
Extended, Blended or Repeating Conjunction Events 


Satellites persistently orbiting in close proximity (e.g., formation flying) can make repeated 
close approaches to one another and produce a probability rate curve R.(f) with multiple blended 
peaks. Figure 5 shows R-(f) curves calculated from an archived OCM involving two close-prox- 
imity satellites. The solid green curve uses all of the 2D P, assumptions and produces a single 
isolated R.(t) peak, as it should. However, the dashed blue curve shows that introducing Keplerian 
two-body motion leads to multiple peaks, caused by repeated encounters between the two satellites. 
Introducing propagated covariances can then change the shape of these peaks. The R-(f) curve es- 
timated using the full Keplerian two-body dynamical model, plotted as the dotted black curve, 
clearly shows multiple, repeating peaks that blend together. Notably, the repeating, nearly periodic 
nature of this probability rate curve is likely an artifact of the relatively simple two-body motion 
model employed here; the effect of orbital perturbations included in more advanced dynamical 
models would likely create more pronounced differences between the sequential peaks. 


Extended, repeating or blended events can be identified by analyzing conjunction durations. 
For instance, all such events contained in this archive analysis can be identified as those with a 
conjunction duration, tA, larger than ~1% of the minimum orbital period for the two objects, Tinin. 
For example, tA/Tinin equals 14% and 65% for the cases shown in Figures 4 and 5, respectively. 
Of the original 80,453 high quality archived conjunctions, 346 have tA/Tinin 2 1%, and produce 
repeating probability rate curves qualitatively similar to that plotted in Figure 5. The remaining 
80,107 OCMs represent single-peaked conjunctions (referred to in this analysis as “isolated” con- 
junctions). 
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Figure 5. A plot of the probability rate for two close-proximity satellites. 
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As mentioned previously, such repeating encounters between close-proximity satellites can po- 
tentially cause Coppola’s formulation to yield inaccurate P. estimates.*'! In these cases, the CARA 
team uses instead P. estimates from a Monte Carlo method, which does not suffer from this limi- 
tation. Notably, in these cases, the P, estimates from Coppola’s method (i.e., those derived by 
performing an integration over the dotted black curve in Figure 5) generally erroneously exceed 
the Monte Carlo P, values. As pointed out by Chan,'! these erroneously large P, values can poten- 
tially even exceed unity, re-emphasizing the limitations of Coppola’s formulation in these situa- 
tions. Despite this, the multi-peaked R.(4) variations still provide a useful indicator of the times 
and frequencies of elevated collision risk for such close-proximity satellites, even though the R. 
integral in Eq. (8) potentially no longer provides an accurate approximation for dP/dt. 


The Frequency of 2D and 3D P, Discrepancies 


After eliminating the extended, repeating or blended conjunctions from the archived data as 
described above, most of the remaining 80,107 cases were found to have negligibly small 3D P. 
values. For instance, the software returns a 3D P, value numerically equal to zero for 44,970 (or 
56%) of these. Only 11,211 (14%) were found to have P. > 10°, and 8,269 (10%) with P. > 10°!°. 
Finally, the CARA team often uses a threshold of P, > 10° for initiating follow-up analyses, and 
only 5,761 (7.2%) of the cases exceeded this threshold. 


Similarly, only 2,674 (or 3.3%) were found to have 3D P, > 10°, which represents the most 
important set for the CARA team. For this set, most 2D and 3D P. estimates were found to be 
relatively close to one another: 1462 (55%) were found to be within 3% of one another, 1,912 
(71%) within 10% of one another, and 2266 (85%) within 30%. Smaller subsets showed larger 
discrepancies: 216 (8.1%) were found to differ by a factor of 2 or more, 149 (5.6%) by a factor of 
3 or more; and 64 (2.4%) by 10 or more. These more discrepant subsets include cases where the 
2D P. estimate is both much greater as well as much smaller than the 3D estimate, which produce 
false alarms and misdetections, respectively.” The misdetections where 3D P. >> 2D P, concern 
the CARA team most. Specifically, the existence of such cases in the archive suggests that a small 
fraction of historical conjunctions with significant risk may have been overlooked because of inac- 
curately low 2D P, estimation. Figures 6 and 7 show plots for two such examples, where the 3D 
P. value exceeds the 2D value by a factor of about four and by several orders of magnitude, respec- 
tively. 


The CARA team continues to analyze the nature and frequency of such large 2D vs. 3D method 
discrepancies for isolated conjunctions. One notable and intriguing feature that such discrepant 
conjunctions seem to show is that the conjunction bounds (to, t1) often do not bracket the time of 
closest approach. This is true for both of the cases plotted in Figures 6 and 7. Other features often 
shown by the most widely discrepant cases are that the peak in the probability rate curve can be- 
come relatively thin and/or occur relatively far from the conjunction midpoint time, both of which 
can be seen in Figure 7. In other discrepant cases not shown here, the (to, t1) bounds do not even 
bracket the R.(t) peak. These effects reinforce the need to implement an adaptive algorithm to cal- 
culate the temporal integration limits and step-sizes used when numerically calculating the integral 
in Eq. (7), as discussed previously. 
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Figure 6. An archived conjunction event where the 2D and 3D Pc values differ by a factor of ~4. 


Other cases of large 2D vs. 3D discrepancies also seem to be characterized by poor quality 
covariance matrices. In this analysis, cases where the relative covariances A(fea) and P(tea) are not 
positive-definite have been eliminated from the analysis. However, some discrepant conjunctions 
that pass this test may still have one (or more) low-quality primary and/or secondary covariances 
among the set {Apj(fea), Pp(tea), As(tca), Ps(tea)}. Again, the CARA team continues to analyze the 
nature and frequency of such discrepancies. 
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Figure 7. An archived event where the 2D and 3D Pc values differ by several orders of magnitude. 


Finally, the specific discrepant cases noted so far have all occurred when using the current ver- 
sion of CARA’s 3D P. software which employs the Keplerian two-body dynamical model and sin- 
gle-Gaussian ECI-state PDF approximations. As stated previously, the current implementation 
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using these approximations only represents a first step, to be followed by a series of upgrades to 
create a more general and accurate implementation of the method. After such upgrades, the specific 
discrepant cases noted here may change or vanish, but likely to be replaced by other similarly dis- 
crepant cases (see Reference 20). However, the plots shown in this initial analysis sufficiently 
illustrate that, for isolated conjunctions, the time of peak collision risk does not necessarily coincide 
with the time of closest approach or even with the conjunction midpoint time, and that the 2D 
approximation can potentially produce very inaccurate results for some conjunctions. 


A Small-P, Screening Test Based on Mahalanobis Distances 


The archive analysis also indicates that a relatively efficient calculation can identify which con- 
junctions will have negligibly small probabilities when using the single-Gaussian PDF approxima- 
tion. This small-P, screening calculation exploits a correlation found to exist between 3D P, values 
and relative-state Mahalanobis distances. More specifically, the P. value for an isolated conjunc- 
tion correlates strongly with the minimum relative position Mahalanobis distance that occurs dur- 
ing the event, as explained below. 


The relative position Mahalanobis distance is a time-dependent quantity defined as follows 
M,(t)=[r? Ate” (11) 


where the time dependence on the relative position vector r and relative position covariance A have 
been suppressed. During an isolated conjunction event, this Mahalanobis distance first decreases 
as the two satellites approach one another, reaches a minimum value, and then increases as they 
separate. The time and value of this minimum Mahalanobis distance can be found using numerical 
methods. Notably, the time that Mp(t) minimizes does not generally coincide with the time of 
closest approach. Analysis of archived conjunctions indicates that it often coincides more closely 
with the conjunction midpoint time, and even more closely yet with the time that the probability 
rate peaks. Analysis also indicates that 3D P, values correlate strongly with the minimum Mp 
value, as shown in Figure 8. The blue crosses show conjunctions with minimum Mp < 10. The red 
crosses show cases with minimum Mp > 10; this larger set comprises about 80% of the conjunc- 
tions, all of which have P. < 3x10°!’, a negligibly small value from the CARA team’s perspective. 


The strong correlation apparent in Figure 8 can be used as a basis for an efficient small-P. 
screening test. Specifically, about 80% of the 3D P, calculations can be circumvented by pre- 
computing the minimum Mahalanobis distance for each conjunction, and eliminating those with 
values larger than 10. The efficiency arises because Mahalanobis distances require significantly 
less computation than 3D P, values. As described previously, the analysis that produced Figure 8 
employs the single-Gaussian PDF and Keplerian two-body motion approximations, but only as a 
first step towards more advanced approaches. These more advanced models will also likely pro- 
duce a correlation like that in Figure 8, which can then be similarly exploited for small-P. screening. 
Furthermore, because these more advanced models will likely require even more intensive compu- 
tation to calculate 3D P. values, performing such screening will potentially become more important 
when analyzing large sets of conjunctions. 
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Figure 8. The correlation between 3D Pc and minimum Mahalanobis distance. 


CONCLUSIONS AND FUTURE WORK 


The implementation of the 3D P, estimation method formulated by Coppola® provides several 
new insights for the conjunction risk assessment community. Many of these insights arise because 
the formulation provides estimates of the probability rate, R-(f). For close-proximity satellites, 
variations in R.(t) can show multiple peaks that repeat or blend with one another, providing a new 
way to view the temporal distribution of risk. For isolated conjunctions, R-(f) analysis provides the 
means to identify and bound times of peak collision risk, which do not always coincide with times 
of close approach or with conjunction midpoint times. Analysis of archived conjunction data 
demonstrates that the commonly used 2D P, approximation can occasionally provide inaccurate 
estimates, including misdetections where the 3D P. value exceeds the 2D value by several orders 
of magnitude. The 3D method also provides the basis for a small-P. screening test that can signif- 
icantly speed risk analysis computations for large numbers of conjunctions. 


As mentioned previously, the CARA team has employed several approximations for the initial 
3D P. software implementation presented here. Future enhancements may entail relaxing the as- 
sumptions behind these approximations and using more advanced algorithms, including: 


1. Adaptive quadrature schemes for the numerical integrations in Eqs. (7) and (8). 
Higher fidelity dynamical models, such as state propagation that includes all relevant orbital 
perturbations. 

3. More advanced covariance propagation such as that derived from an Equinoctial orbital 
element description.!?° 


4. More realistic state PDF representations, such as the GMM expansion from Eq. (2).!2!415 


By implementing these upgrades, we believe that the 3D P. method has the potential to provide 
substantially more realistic and accurate conjunction risk assessments. 
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